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Chapter 1 

Nonequilibrium Stationary Solutions of Thermostated 
Boltzmann Equation in a Field 

F. BonettcQ, J.L. LebowitsQ 



Dedicated to Leopoldo Garcia Colin on the occasion of his eightieth 

birthday. 



We consider a system of particles subjected to a uniform external force 
E and undergoing random collisions with "virtual" fixed obstacles, as 
in the Drude model of conductivity [1]. The system is maintained in a 
nonequilibrium stationary state by a Gaussian thermostat. In a suitable 
limit the system is described by a self consistent Boltzmann equation for 
the one particle distribution function /. We find that after a long time 
/(v, t) approaches a stationary velocity distribution /(v) which vanishes 
for large speeds, i.e. /(v) = for |v| > v max (E), with w ma:E (E) ~ |E| _1 
as |E| — > 0. In that limit /(v) ~ exp(— c|v| 3 ) for fixed v, where c 
depends on mean free path of the particle, /(v) is computed explicitly 
in one dimension. 
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1.1. Introduction 

Nonequilibrium stationary states (NESS) in real systems must be main- 
tained by interaction with (effectively infinite) external reservoirs. Since a 
complete microscopic description of such reservoirs is generally not feasible 
it is necessary to represent them by some type of modeling [2] . 

However, unlike systems in equilibrium, which maintain themselves 
without external inputs and for which one can prove (when not inside a 
coexistence region of the phase diagram) that bulk behavior is independent 
of the nature of the boundary interactions, we do not know how different 
microscopic modeling of external inputs, affects the resulting NESS. 

One particular type of modeling dynamics leading to NESS is via Gaus- 
sian thermostats, see [3]. Analytical results as well as computer simulations 
have shown that the stationary states produced by these dynamics behave 
in many cases in accord with those obtained from more realistic models [4] . 
This has led us to continue our study of the NESS in current carrying ther- 
mostated systems. In our previous work, see [3, 5] we carried out extensive 
numerical and analytical investigations of the dependence of the current 
on the electric field for a model system consisting of N particles with unit 
mass, moving among a fixed periodic array of discs in a two dimensional 
square A with periodic boundary conditions, see Fig. 1. They are acted on 
by an external (electric) field E parallel to the x-axis and by a "Gaussian 
thermostat". (The discs are located so that there is a finite horizon, i.e. 
there is a maximum distance a particle can move before hitting a disc or 
obstacle). 

The equations of motion describing the time evolution of the positions 
qi and velocities Vj, i = 1,...,N, are: 

q» = v 2 q; = {q ilX ,qi,y) G A' (1.1) 
v; = E-a(J,U)v l +F obs (q l ) (1.2) 

where 

T F 1 N 1 N 

«(j.to = ^. J = ^E v <> u = jf^ (L3) 

i=l i=l 

Here A' = A\V, with V the region occupied by the discs (obstacles) and 
F b s represents the elastic scattering which takes place at the surface of 
the obstacles. The purpose of the Gaussian thermostat, represented by 
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Fig. 1.1. General billiard structure with three particles shown. 

the term a(J, U)v in Eq. (jl.ip . is to maintain the total kinetic energy 
1/2 YliLi v ? constant, i.e. U = v^. It also has the effect of making the flow 
$t generated by Eq. (jl.ip on the (AN — 1) dimensional energy surface non 
Hamiltonian when E ^ 0. Another effect of the thermostat is to effectively 
couple all the particles in a mean field way a(J,U), depending only on 
the total momentum of the particles. Note that this is the only coupling 
between the particles in this system. 

Our main interest is in the NESS of this model system. To get some 
analytical handle on the form of the NESS we also investigated numerically 
a model system in which the deterministic collisions with the obstacles 
are replaced by a stochastic process in which particle velocities get their 
orientations changed at random times, independently for each particle [5]. 
This yields, in the limit N — > oo, a self consistent Boltzmann equation. In 
the spatially uniform case, the only one we shall consider here, the equation 
takes the form 

^/(v, t) + ^ [(E - M v)/(v, t)] = -Xj (/(v, t) - /(v', t)) dn 

(1.4) 

where v' = v — 2n(n • v), n is a unit vector and [i is determined self 
consistently by requiring that 
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d_ 

dt 



y^/(v,i)dv = -E.<v) +A1 (v 2 )=0 (1.5) 



or 



H= — X- (1.6) 



E-(v) 

(v 2 ) 

This ensures that (v 2 ) is independent of time, i.e. J v 2 f(v,t)dv = 
J v 2 /(v,0)<iv. When E = also /i = and /(v, t) will approach, as 
t — ► oo, the function /(|v|,0) = / /(v, 0)cm, where 5(1) is the area of 
the unit sphere in n dimensions. There will thus not be a unique stationary 
solution of Eq. (|1.4[) . The situation is different when E ^ 0. In that case 
limt^oo /(v, t) = /(v) independent of /(v, 0). It is this NESS of Eq. (fl~4|) 
which we shall study here. 

The outline of the rest of the paper is as follows. In section 2 we show 
that the stationary solution vanishes for |v| > |E|//i. In section 3 we 
obtain an exact expression for the steady state solution of this equation 
in one dimension. In section 4 we investigate the small field limit of this 
stationary distribution. Finally in section 5 we find the small field limit of 
the distribution in two dimensions. 

1.2. Properties of Eq. (I1.4[) in arbitrary dimension 

Here we investigate some general properties of Eq. (|1.4j) in the steady state. 
Observe that the right hand side of the equation preserves v 2 . On the other 
hand we have 



E 

(E - (mv) • v < if Ivl > -! — L (1.7) 

A 4 

Thus if we define 



we get that 



d 
dt 



F(V,t)= f f(y,t)dv 

J\v\<V 

F(V, t) = - f ([(E • v) - /iy 2 ]/(v, <)dv < 

J\v\ = V 



(1.8) 



(1.9) 

<-S(V)(\E\-ti\V\) sup /(v,i) (1.10) 

|v|=V 
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where S(V) is the surface of the sphere of radius V. Since in the steady 
state we must have F(V, t) = it follows that /(v) = if v > |E|//x, where 
fjj is now the stationary self consistent value. When |E| — > 0, /i will decrease 
as | E | 2 , as long as A > 0, and so the maximum speed will grow as | E | 1 . 

Observe that, as noted before, the limit for E — > is very different from 
that obtained when E = and \i = 0. This behavior is due to the fact that 
the collision term preserves |v| 2 . In particular it would fail if one also adds 
a non linear Boltzmann particle-particle collisions term to the right hand 
side of Eq. I|1.4p . In that case the stationary distribution would be positive 
for all v. We expect however that Eq. (| 1 .4|) would describe the behavior of 
the system for a long time when the effective rate of inter-particle collisions 
is small compared to the rate of collisions with the obstacles. 

1.3. The Boltzmann Equation in one dimension 

We consider the one dimensional version of Eq. (|1.4p . This takes the form: 

§- t f(v, t) + ^ [{E - H/M)] = -AM (f(v, t) - /(-«,*)) (1.11) 

with /j, determined self consistently by the requirement that the average 
energy remain fixed, see Eq. (|1.6p . We look for a stationary solution of 
Eq. (|1.4p in the form 

f(v)=4>(v) + Eil>(v) (1.12) 

with 

<p{v) = (t>{-v) 
tP(- v ) = -i/)(-v) (1.13) 

and 4> and ip defined for v > 0. The time independent Eq. (jl.lip takes the 
form 



^ \fiv<l>(v) - E 2 i,(v)] = 
V d 

— Wv) - fivip(v)] = -2\vMv) (1.14) 
ov 

where /x now has its steady states value given by Eq. (|1.6[) . Using the fact 
that iji is odd the first equation implies that 
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iftv) = g^(v). (1.15) 

Setting /i = vE 2 we have tp(v) — vv4>{v) and the second equation in 
Eq. (II. 14p becomes 

^- [(1 - v 2 E 2 v 2 )<f>{v)} = -2\vv 2 <f>(v) (1.16) 

Observe that since f(v) > we have from the definition, Eq. (| 1 . 1 3|) . that 
4>{v) > 0. It follows from Eq. (|1.16|) that, if <f>(v) > for v > then we 
will have lim^oo <fi — oo, so we must have 

cf>(v) =0 for v > — (1.17) 

vE 

This is consistent with the general result in the previous section. Setting 
$(v) = (1 - v 2 E 2 v 2 )(t>{v) (1.18) 

Eq. (|1.16p becomes 

d , $» 2\vv 2 2A 2A 1 , . 

— loe $(v) = — — = = (1 19) 

dv S ( ' $(v) 1 - v 2 E 2 v 2 vE 2 vE 2 1 - v 2 E 2 v 2 K ' ' 

This can be solved to give 
or 

^ = T^W (iT^J e ^ (L21) 

for v < and otherwise. We finally get, using Eq. (| 1 . 1 3[) . that for «£l, 



l - V 1 + 

for < and otherwise. C is a normalization constant. 

Observe that at v = (vE)^ 1 the solution is singular if A < v 2 E 3 . This 
means that the field is strong enough to push the velocity of the particle 
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close to its limiting value before a collision happerQ. Furthermore, when 
the field goes to we do not obtain a Maxwellian distribution. To be 
sure, there is no reason why this should happen since we do not have any 
mechanism (such as inter-particle collisions) which would tend to bring the 
system to equilibrium, as one would expect to be the case in a realistic 
physical model. 

It follows from Eq. (fTTTC|) that 

fv><p(v)dv = £- v (1.25) 

so that the average kinetic energy in the steady state is (v 2 ) = C j\v = Vq 
and the average current is (v) = CE/Xv. 



1.4. Small field limit 

Equation (|1.22[) can be rewritten in the form: 



f(v) 



C 



1 - vEv 



exp 



A 



v 2 E 3 



(log (1 - uE\v\) - log (1 + uE\v\)) 



(1.26) 

This expression is clearly analytic in E for E < l/u\v\. In particular we can 
compute the zero order term by expanding log (1 — vE\v\) — log (1 + vE\v\) 
in term of Clearly only terms odd in E will appear. The term linear 



in E cancel out the exponent 



2A 



The term proportional to E gives, 



once multiplied by X/v 2 E 3 , a term 2/3A^|w| 3 . This gives us: 



f(v) = C(l + Evv) exp f~\ vv 3 j + o{E) (1.27) 



a Note that if we let A — ► than the singularity becomes non integrable and: 

lhn^) = i 5 ( H --L) (1.23) 

so that 

f(v)=S (v--L^ (1.24) 

with (v 2 ) = (1/uE) 2 = Vq specified a priori, so that u = 1/Evq or = E/vo. We 
also have (v) = vq. This is exactly what happens in the TV-particles d-dimcnsional 
thermostated case without collisions. All particles end up moving parallel to the field 
with the same speed independent of the strength of |E| > 0. 
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where the term linear in E comes from the prefactor 1/(1 — vEv). Here v 
is the limiting value of v which, given the (initial) kinetic energy Vq/2, is 
given by 

v=-^r- + 0{E 2 ) (1.28) 

with 

K = - 3^ (1.29) 

8 7T2 

We can also obtain the behavior of f(v) when E tends to by setting 
E = in Eq. (jTTBJ). This gives: 

— (j>{v) = -2\vv 2 (t>{v) (1.30) 
whose solution is found to be 

4>{v) = Ce"* W (1.31) 
This agrees with the result in Eq. (|1.27[) when E — ► 0. 

1.5. The general case for small E 

In the two dimensional case the stationary version of Eq. (jl.4p takes the 
form 

A [(E - HE| 2 v)/(v)] = -\ I ^ (/(v) - /(v')) rfn (1.32) 

where v' = v — 2n(n • v), n is a unit vector and we have set fj, — ^|E| 2 . 

We do not have an exact solution of this equation for arbitrary field 
strength. We can however still find the small |E| behavior of /(v) as we 
did before for the one dimensional case. To do this we first expand /(v) in 
harmonics of the angular part of v. This means that that we write: 



/( v ) = ]T cos(n0) (1.33) 

n=0 

where v = (v cos(0), v sin(0)). We have taken the field to be in the re- 
direction, E = {E, 0), so that no term in sin(#) appears due to the symmetry 
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of the system with respect to the x axis. Substituting in Eq. (|1.32p we 
obtain: 

- E 2 vvd v f n {v) - 2E 2 uf n (v) + (1.34) 

rpx /n-i(u) + fn+i(v) „{n- l)/„-i - l)/„+l 
+ Ed v E = 



An 2 - 1 



1 \vf n iv) 



where we have set /_i = 0. Since we expect that the distribution /(v) 
will depend only on |v| when E — > we can assume that f n (v) = Eg n (v) 
for n ^ 0. When n = 0, the right hand side of Eq. (|1.34[) vanishes. The 
equation thus becomes, after symplifying a common E 2 factor, 

vvd v fd{v) + 2vf (v) - \d v g x (v) - ^-gi{v) = (1.35) 
while the leading term in E of the equation for n = 1 is 

d v f (v) = ~\v 9l (v). (1.36) 
So that, at first order in E, we have: 

/(v) = C(l + 2^E ■ v) exp ^-^|v| 3 ) (1.37) 

This is the same form as Eq. (|1.27[) with the limiting value for v given by 
Eq. (OH) with 

13^ 

9r(l) 3 

1.6. Concluding remarks 

A similar analysis shows that the small |E| behavior of the NESS described 
by Eq. (|1.4[) will be of the same form as that given by Eq. (|1.37p in all 
dimension d, in particular the zero order term will be an exponential in 
— |v| 3 . We do not however have any intuitive heuristic explanation for this 
behavior on the more microscopic many particles level. 

Let us go back now to our original problem, the NESS for the N particles 
system in an electric field and Gaussian thermostat, described by Eqs. 
(|1.1[) - (|1.3[) . As mentioned in the introduction, we expect that in the limit 
N — > og this NESS will agree with our results for the stationary distribution 
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of the Boltzmann equation (|1.4[) i.e. there will be a cut-off in the particle 
speed independent of N when N becomes very large. This indicates that 
the NESS for the original model will be, for large N, concentrated on those 
parts of the dN — 1 dimensional sphere in which the particle speeds are 
closer to each other than they would for a uniform distribution on the 
sphere. The latter is of course what gives rise to a Maxwellian distribution 
with the same mean energy. This is consistent with what we found in [5], 
see in particular section III there. 
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